

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z analyst specif


%% estimated mean utilities
UStar = payoffs(thetaStar,X,Z,YY,reb,mpow);


%% counterfactuals
load('counterfactualRU_replication.mat','srng_cfactRU')
rng(srng_cfactUS)
NS = 2e+3;
YY4 = YY(YY(:,6)==0,1:5); % 4-mpow cw-intervention strategy profiles
timeRU=0;

% RU never intervenes
VnoRU_R = UStar(YY(:,5)==0,1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VnoRU_US = UStar(YY(:,5)==0,2,1:N) + randn(size(YY4,1),NS,N);
VnoRU_UK = UStar(YY(:,5)==0,3,1:N) + randn(size(YY4,1),NS,N);
VnoRU_France = UStar(YY(:,5)==0,4,1:N) + randn(size(YY4,1),NS,N);
VnoRU_Russia = UStar(YY(:,5)==0,5,1:N) + randn(size(YY4,1),NS,N);
VnoRU_China = UStar(YY(:,5)==0,6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
noRUintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameRU',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "UK" "France" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VnoRU_US(2:end,b,n),VnoRU_UK(2:end,b,n),VnoRU_France(2:end,b,n),VnoRU_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_UK=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_France=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_UK.^YY4(2:end,3)).*((1-sigma_UK).^(1-YY4(2:end,3))).*...
                (sigma_France.^YY4(2:end,4)).*((1-sigma_France).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VnoRU_R(2:end,b,n)'*EQ(:,m);
        end
        noRUintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeRU=timeRU+toc;
end
[noRUDMES,noRUintervEq,noRUeqPayoffs] = createDMES(noRUintervEq,VnoRU_R,VnoRU_US,VnoRU_UK,VnoRU_France,VnoRU_Russia,VnoRU_China,[YY4(:,1:4),zeros(size(YY4,1),1),YY4(:,5)],specif);

% RU always intervenes
VallRU_R = UStar([1;find(YY(:,5)==1)],1,1:N) + randn(size(YY4,1),NS,N); % entry payoffs of rebel group
VallRU_US = UStar([1;find(YY(:,5)==1)],2,1:N) + randn(size(YY4,1),NS,N);
VallRU_UK = UStar([1;find(YY(:,5)==1)],3,1:N) + randn(size(YY4,1),NS,N);
VallRU_France = UStar([1;find(YY(:,5)==1)],4,1:N) + randn(size(YY4,1),NS,N);
VallRU_Russia = UStar([1;find(YY(:,5)==1)],5,1:N) + randn(size(YY4,1),NS,N);
VallRU_China = UStar([1;find(YY(:,5)==1)],6,1:N) + randn(size(YY4,1),NS,N);
filen=1:NS;
allRUintervEq=cell(NS,N);
for n=1:N
    tic
    parfor b=1:NS
        fn=filen(b);
        fname=['gameRU',num2str(fn),'.nfg'];
        fileID=fopen(fname,'w');
        fprintf(fileID,'NFG 1 R "Major-Powers Intervention"\n',1); %#ok<CTPCT>
        fprintf(fileID,'{ "US" "UK" "France" "China" } { 2 2 2 2 }\n',1); %#ok<CTPCT>
        formatSpec='%12g ';
        fprintf(fileID,formatSpec,reshape([VallRU_US(2:end,b,n),VallRU_UK(2:end,b,n),VallRU_France(2:end,b,n),VallRU_China(2:end,b,n)]',1,4*2^4)); %#ok<*PFBNS>
        fclose(fileID);        
        [~,NE]=system(['gambit-enumpoly ',fname]);
        findNE=strfind(NE,'NE');
        EQ=zeros(2^4,length(findNE));
        EV_R=zeros(1,length(findNE));
        for m=1:length(findNE)
        	sigma_US=str2double(NE(findNE(m)+12:findNE(m)+19));
            sigma_UK=str2double(NE(findNE(m)+30:findNE(m)+37));
            sigma_France=str2double(NE(findNE(m)+48:findNE(m)+55));
            sigma_China=str2double(NE(findNE(m)+66:findNE(m)+73));
            EQ(:,m)=(sigma_US.^YY4(2:end,2)).*((1-sigma_US).^(1-YY4(2:end,2))).*...
                (sigma_UK.^YY4(2:end,3)).*((1-sigma_UK).^(1-YY4(2:end,3))).*...
                (sigma_France.^YY4(2:end,4)).*((1-sigma_France).^(1-YY4(2:end,4))).*...
                (sigma_China.^YY4(2:end,5)).*((1-sigma_China).^(1-YY4(2:end,5)));
            EV_R(m)=VallRU_R(2:end,b,n)'*EQ(:,m);
        end
        allRUintervEq{b,n}=[EV_R;EQ];
    end
	toc
    timeRU=timeRU+toc;
end
[allRUDMES,allRUintervEq,allRUeqPayoffs] = createDMES(allRUintervEq,VallRU_R,VallRU_US,VallRU_UK,VallRU_France,VallRU_Russia,VallRU_China,[YY4(:,1:4),ones(size(YY4,1),1),YY4(:,5)],specif);


%%
clearvars -except N NS srng_cfactRU YY4 timeRU VnoRU_R VnoRU_US VnoRU_UK VnoRU_France VnoRU_Russia VnoRU_China ...
     VallRU_R VallRU_US VallRU_UK VallRU_France VallRU_Russia VallRU_China noRUintervEq noRUDMES noRUeqPayoffs allRUintervEq allRUDMES allRUeqPayoffs
save('counterfactualRU.mat')



